%function [h,pValue,stat,cValue] = chidtest(sP,eP,sN,eN,dof)
%CHIDTEST chi-square distance test of two matrices.

% It is to test whether the transition matrices estimated from each of the
% T subsamples differ sigificantly from the matrix estimated from the
% entire sample.

% Input Arguments:
%
%   sP - subsample transition probability matrix. (txixj) dimention. Initial rating i = 1, 2, ..., d-1 and
%   the end rating j = 1, 2, ..., d (including default)
%
%   eP - entiresample transition probability matrix. Initial rating i = 1, 2, ..., d-1 and
%   the end rating j = 1, 2, ..., d (including default)
%
%   sN - subsample transition number matrix. (txixj) dimention. Initial rating i = 1, 2, ..., d-1 and
%   the end rating j = 1, 2, ..., d (including default)
%
%   eN - entiresample transition number matrix. Initial rating i = 1, 2, ..., d-1 and
%   the end rating j = 1, 2, ..., d (including default)
%
%   dof - Degree-of-freedom parameters for the asymptotic chi-square
%       distributions of the test statistics. Elements of dof are positive
%       integers equal to (d-1)*(T-1). If dof is a scalar, it is expanded to a
%       vector with length equal to the number of tests. If dof is a
%       vector, it must have length equal to the number of tests.

% Output Arguments:
%
%   h - Vector of Boolean decisions for the tests, with length equal to the
%       number of tests. Values of h equal to 1 indicate rejection of the
%       null, restricted model in favor of the alternative, unrestricted
%       model. Values of h equal to 0 indicate a failure to reject the
%       restricted model.
%
%   pValue - Vector of p-values of the test statistics, with length equal
%       to the number of tests.
%
%   stat - Vector of test statistics, with length equal to the number of
%       tests.
%
%   cValue - Vector of critical values for the tests, determined by alpha,
%       with length equal to the number of tests.

% Load data
clc;
clear;
load TP06_Y_RJ_fin_R5_G250_1000k_MR7.mat %final_result=zeros(13,10,13);
load Prob_std_TP06YRJfinR5G2501000kMR7.mat % prob_std=zeros(13,10,12);
load Avg_Prob_std_TP06YRJfinR5G2501000kMR7.mat %stdprob=zeros(1,9,10);  
load Counts_TP06YRJfinR5G2501000kMR7.mat; %count=zeros(1,10,13); total counts

sP = zeros(12, 7, 8); % year 1998 -2009, total 12yr; rating 1.5-4.5 and default total 8 categories.
sN = zeros(12, 7, 8);
eP = zeros(7, 8);
filename='TP06_Y_RJ_fin_R5_G250_1000k_MR7_ChiSqureTest.xls';




% Collect data
sP(:, :, 1:7) = prob_std(2:13, 3:9, 3:9);
sP(:, :, 8) = prob_std(2:13, 3:9, 11);
sN(:, :, 1:7) = final_result(2:13, 3:9, 3:9);
sN(:, :, 8) = final_result(2:13, 3:9, 11);
eP(:, 1:7) = stdprob(1, 2:8, 2:8);
eP(:, 8) = stdprob(1, 2:8, 10);


% Check alpha input set default:
alpha = 0.05;
stat = zeros(8,1); % the first 7 is for different rating, and last one is for the sum of time.
dof = zeros(8,1);
pValue = zeros(8,1);
h = zeros(8,1);
cValue = zeros(8,1);

num = 0;

% Perform the tests:
for i = 1:7
for t=1:12
    for j=1:8
        stat(i)= sum(sN(t,i,1:8)).*(sP(t,i,j) - eP(i,j)).^2./eP(i,j) + stat(i);
        if (sum(sN(t,i,1:8))* eP(i,j))< 5
            num = num +1;
        end
    end
end
dof(i) = (size(sP,1)-1)*(size(sP,3)-1);
pValue(i) = 1-chi2cdf(stat(i),dof(i));
h(i)= (pValue(i) <= alpha);
cValue(i) = chi2inv(1-alpha,dof(i));
stat(8) = stat(i) + stat(8);
end

dof(8) = (size(sP,1)-1)*(size(sP,3)-1)^2;
pValue(8) = 1-chi2cdf(stat(8),dof(8));
h(8) = (pValue(8)<= alpha);
cValue(8) = chi2inv(1-alpha,dof(8));

sheetlist=strcat('Parameters');
d = {'h','pValue','stat','cValue','dof','Num of less than 5'; h, pValue, stat, cValue, dof, num};
xlswrite(filename, d, sheetlist)



% for t = 1 : 12
%     for i = 1 : 7
%         for j = 1 : 8
%         stat = sum(sN(t,i,1:8)).*(sP(t,i,j) - eP(i,j)).^2./eP(i,j) + stat;
%         if (sum(sN(t,i,1:8))* eP(i,j))< 5
%             num = num +1;
%         end
%         end 
%     end
% end 

%     
% 
% dof = (size(sP,1)-1)*(size(sP,3)-1)^2;
% pValue = 1-chi2cdf(stat,dof);
% h = (pValue <= alpha);
% cValue = chi2inv(1-alpha,dof);



